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Abstract 

In a recent paper Q a new powerful method to calculate Feynman diagrams 
was proposed. It consists in setting up a Taylor series expansion in the 
external momenta squared. The Taylor coefficients are obtained from the 
original diagram by differentiation and putting the external momenta equal 
to zero. It was demonstrated that by a certain conformal mapping and 
subsequent resummation by means of Pade approximants it is possible to 
obtain high precision numerical values of the Feynman integrals in the 
whole cut plane. The real problem in this approach is the calculation of 
the Taylor coefficients for the arbitrary mass case. Since their analytic 
evaluation by means of CA packages uses enormous CPU and yields very 
lengthy expressions, we develop an algorithm with the aim to set up a 
FORTRAN package for their numerical evaluation. This development is 
guided by the possibilities offered by the formulae manipulating language 
FORM §. 



* to be published in the proceedings of the workshop on Computer Algebra in Science and Engineering, 
held at Zcntrum fiir Interdisziplinare Forschung (ZiF), Bielefeld, Germany, August 1994; World Scientific, 
J. Fleischer, J.Grabmcicr, F.W. Hehl and W.Kuchlin, editors. 

t Supported by Bundcsministerium fiir Forschung und Technologic. 

t On leave of absence from Joint Institute for Nuclear Research, 141980 Dubna, Moscow Region, 
Russian Federation. 



1. Introduction 

Standard-model radiative corrections of high accuracy have obtained growing at- 
tention lately in order to cope with the increasing precision of LEP experiments Q. In 
particular two-loop calculations with nonzero masses became relevant Q. While in the 
one-loop approach there exists a systematic way of performing these calculations ||, in 
the two-loop case there does not exist such a developed technology and only a series of 
partial results were obtained ||, but no systematic approach was formulated. 

Our approach consists essentially in performing a Taylor series expansion in terms 
of external momenta squared and analytic continuation into the whole region of kinemati- 
cal interest. Simple as this may sound, there are some unexpected methodical advantages 
compared to other procedures. 

Considering a Taylor series expansion in terms of one external momentum squared, 
q 2 say, the differential operator by the repeated application of which the Taylor coefficients 
are obtained, subsequently setting q — 0, is 

d 2 

Dg= dq~dq^- (1) 

Such expansions were considered in |5J, Pade approximants were introduced in || and 
in Ref. [Q] it was demonstrated that this approach can be used to calculate Feynman 
diagrams on their cut which, concerning physics, is the most interesting case. The above 
Taylor coefficients are essentially "bubble diagrams", i.e. diagrams with external mo- 
menta equal zero. They are essentially the same (after partial fraction decomposition) 
for two-point, three-point, . . . functions for a given number of loops and we stress that 
it is indeed a great technical simplification to have to perform integrals only for exter- 
nal momenta equal zero even if these integrals contain now arbitrary high powers of the 
scalar propagators. For their calculation recurrence relations are quite effective ( || , [|HJ] , 
see Sect. 6). On this basis we develop our algorithm. 

2. Expansion of three-point functions in terms of external momenta squared 

Here we have two independent external momenta in d = 4 — 2e dimensions. The 
general expansion of (any loop) scalar 3-point function with its momentum space repre- 
sentation C(pi,p2) can be written as 

oo oo 

C(pi,ft)= E a lmn { V \)\vlT(VlV2Y = E E *lmn(p 2 l) l (pt) m (PlP2) n , (2) 
(,m,n=0 L=0 l-\-m+n=L 



where the coefficients ai mn are to be determined from the given diagram. They are 
obtained by applying the differential operators = q^-^Jf several times to both sides 

of @. 

This procedure results in a system of linear equations for the a/ mn . For fixed L 
(see equation (||)) we obtain a system of (L + 1)(L + 2)/2 equations of which, however, 
maximally [L/2] + 1 couple ([x] standing here for the largest integer < x). These linear 



equations are easily solved with REDUCE [11], e.g., for arbitrary d. 



For the purpose of demonstrating the method, we confine ourselves to the case 
Pi = p\ = 0i which is e.g. physically realized in the case of the Higgs decay into two 
photons (H — > 77) with p\ and P2 the momenta of the photons. In this case only the 
coefficients aoon are needed. They are each obtained from a "maximally coupled" system 
of [n/2] + 1 linear equations. Solving these systems of equations we obtain a sequence of 
differential operators (-D/'s) which project out from the r.h.s. of (Q) the coefficients aoo n : 

Df - [n V +1 (-^-T^ + n-Qrtrf-l) 

Dfo0n ~ ti 2r(*)r(n - 2, + 3)r(n + d - 2)r(n + d/2) ( 12) ( 11 22) ' (3) 
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where the sum of the exponents of the various □' s is equal n. Applying the operator -D/oon 
to the (scalar) momentum space integral C(pi,p2) and putting the external momenta equal 
to zero, yields the expansion coefficients aoon- 

In the two-loop case we consider the scalar integral (k% = k% — k 2 , see also Fig. 1) 



C(mi, • • ■ ,m e ;p 1 ,p 2 ) 

72)2 J ' 



d A k 1 d i k 2 



(4) 
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Figure 1: Planar and non-planar scalar vertex diagrams and their kinematics 

Introducing the abbreviations C\ = k\—m\, c 2 = k\—m\, C3 = k\— wij, C4 = k\—m\ 
and C5 = k\ — mf, c§ = k\— wig, we have (C5 and c§ do not enter the differentiation since 
for the planar as well as for the non-planar diagram they occur in @) as such) 

2' 



(iTT 2 ) 2 a on 



- [ d%d 4 k 2 F n - 

1 J Ci Co Ca C 



n + 1 J c\ c 2 c 3 C4 c 5 c 6 

where for the planar diagram (k% = k 2 in C4 and k% = k\ — k 2 in cq) 

n n 



u=0 



u'=0 



and 



fea) = E aZiklT-^'^iklYih k 2 y + »'~ 2 ^ 

0<2fi<v+v'<n+n 



a^, being rational numbers with the properties 



nfi nfi 

a vv' ~ a v'v 



and EX£ = 1 > 



(5) 



(6) 
(7) 



The above is now essentially the basis for the algorithm we are going to develop: 
to calculate the Taylor coefficients of integrals like (|) according to (|J) - (|8|) and to reduce 
them to integrals of the type 
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^(a,/3, 7 ,m 1 ,m 2 ,m 3 ) = (-l)(°+/^> / — — — * k fJ%. (9) 

J (k( — mfp(k 2 — m 2 p{{ki — k 2 ) 1 — m^p 

which in turn will be reduced by means of recurrence relations to Vb(1, 1, 1, mi, 1712, WI3). 
Accordingly our algorithm is performed in the following three steps: 

• First of all the coefficients a^, in (0) and the corresponding ones for the non-planar 
diagram can be evaluated in terms of multiple sums over T-functions. While (^) 
has been obtained in by inspection of FORM output and some lower coefficients 
could be read off explicitly, in Sect. 4 we give a proof of this representation by 
construction, which also yields the 



From (|9]) it becomes clear that the numerator scalar products in (|7D must be elim- 
inated and/or a partial fraction decomposition of products of scalar propagators 
with the same integration momentum fcj, k 2 or k% but different masses must be 
performed. Substituting, e.g., k\ = c$ + mf (i=l,2) one cancels k\. Similarly one 
proceeds for k 2 and k\. For k\k 2 in ([^D one writes for the planar diagram 

hk 2 = -{k\ + k\ - ml) - lc 6 = h 2 - ^c 6 (10) 
and by stepwise reducing higher powers ( v + v' — 2ji = A ) 

2\k 1 k 2 ) x = (k 2 ) x - Uk 2 ) x - X + 2k l k 2 {k 2 ) x ~ 2 + ... + (2k 1 k 2 ) x ~ 1 } c 6 . (11) 



In the second term cq cancels after insertion into (||) so that only factorized one-loop 
contributions are obtained from it (i.e. the integral @ factorizes into integrals over 
ki and k 2 separately). Moreover, in the square bracket of (0) only even powers of 
k\k 2 contribute after integration. The "genuine" two- loop contributions are then 
obtained by replacing k\k 2 in (^) by \{k\ + k\ — mj?) according to the first term 
in (|TTD. The problem of cancelling k\k 2 is more complicated for the non-planar 



diagram, as will be demonstrated in Sect. 5. 

The evaluation of the integrals @ is supposed to be performed in terms of recurrence 
relations, thus reducing them to known two-loop "master" integrals, as will be 
described in Sect. 6. If one of the indices a,/3, 7 in (|) is < (i.e. in our notation 
the corresponding scalar propagator occurs with positive power in the numerator), 
the integral can be expressed again in terms of factorized one-loop integrals and 
simple explicit representations can be found in this case. An example is also given 
at the end of Sect. 6. 



3. The method of analytic continuation 

Before entering the details of the calculation, we want to give a motivation for the 
above algorithm, i.e. the main interest in Feynman diagrams is for the values on their cut 
and we have to demonstrate how to obtain these from the Taylor expansion. 

Assume, the following Taylor expansion of a scalar diagram or a particular ampli- 
tude is given: 

00 

C( Pl , P2 ,...) = Y.*mv m = f(y) (12) 

m=0 
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and the function on the r.h.s. has a cut for y > yo- In the above case of H — > 77 one 

introduces y = t 2 -? with g 2 = (pj — p 2 ) 2 as adequate variable with y = 1. 

Our proposal for the evaluation of the original series is in a first step a conformal 
mapping of the cut plane into the unit circle and secondly the reexpansion of the function 
under consideration into a power series w.r.to the new conformal variable. A variable 
often used lOI is 



— v**. (13) 

V yo 

Considering it as conformal transformation, the y-plane, cut from y to +00, is 
mapped into the unit circle and the cut itself is mapped on its boundary, the upper semi- 
circle corresponding to the upper side of the cut. The origin goes into the point uj = 0. 

After conformal transformation it is suggestive to improve the convergence of the 
new series w.r.to to by applying one of the numerous summation methods ||13||,||14| most 
suitable for our problem. We obtained the best results with the Pade method and partially 
also with the Levin v transformation. The expansion of f(y) in terms of uj is: 

00 

/(y(w)) = 5>V., (14) 

s=0 



where 



n=l 



r(2n)r(s-n + l)' 



Eq.([14|) will be used for the analytic continuation of / into the region of analyticity 
(y < yo; observe that the series ( |P2"D converges for \y\ < y only) and in particular for the 
continuation on the cut (y > y ). In this latter case we write 

uj = exp[z£ (?/)], with cos£ = — 1 + 2— (16) 

and hence 

f{y) = ao + J2 0n exp irt£(y) (17) 

n=l 

In any case we have \u\ < 1 and we will show in the following how to sum the 
above series. 

Pade approximations are indeed particularly well suited for the summation of the 
series under consideration. In the case of two-point functions they could be shown in 
several cases (see e.g. ) to be of Stieltjes type (i.e. the spectral density is positive). 
Under this condition the Pade's of the original series (|1^) are guaranteed to converge 
in the region of analyticity. For the three-point function under consideration (H — ► 
77), however, the obtained result shows that the series is not of Stieltjes type (i.e. the 
imaginary part changes sign on the cut). 

Having performed the above uj transformation (0), however, it is rather the Baker- 
Gammel- Wills conjecture ( see fl5| ), which applies. 

A convenient technique for the evaluation of Pade approximants is the e-algorithm 
of |T3|. In general, given a sequence {S n \n = 0, 1,2, . . .}, one constructs a table of ap- 



proximants using 

T(m,n) = T(m- 2,n + l) + l/{T(m - l,n+ 1) - T(m - l,n)}, (U 
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with T(0,n) = S n and T(—l,n) = 0. If the sequence {S n } is obtained by successive 
truncation of a Taylor series, the approximant T(2k,j) is identical to the [k + j/k] Pade 
approximant derived from the first 2k + j + 1 terms in the Taylor series. 

We present results for the two-loop three-point scalar ( planar) integral with the 
kinematics of the decay H — > 77. We study the integral (|]) with m 6 = and all other 
masses m; = m t (i = 1, ..,5). In this special case all Taylor coefficients can be expressed 
in terms of T -functions. For a list of the first coefficients aoo n (= a n ', n = 0, 28) see [p]]. 



Table 1: Results on the cut (q 2 > Am 2 ) in comparison with |TB| 





[10/10] 
Re 


Im 


[14/14] 
Re 


Im 


Ref.Q 
Re 


Im 


4.01 


11.926 


12.66 


11.935 


12.699 


11.9347(1) 


12.69675(8) 


4.05 


5.195 


10.48 


5.1952 


10.484 


5.1952(1) 


10.4836(4) 


4.10 


2.6624 


9.095 


2.66245 


9.0955 


2.66246(2) 


9.0954(2) 


4.20 


0.5161 


7.4017 


0.516039 


7.401640 


0.51604(5) 


7.40163(4) 


4.50 


- 1.42315 


4.77651 


- 1.42315097 


4.77651003 


- 1.423122(9) 


4.776497(9) 


5.0 


- 1.985805 


2.758626 


- 1.985804823 


2.758626375 


- 1.98580(2) 


2.758625(2) 


6.0 


- 1.7740540 


1.1232494 


- 1.774053979 


1.123249363 


- 1.77405(1) 


1.123250(6) 


7.0 


- 1.4192404 


0.4807938 


- 1.419240377 


0.4807938045 


- 1.419240(5) 


0.480794(9) 


8.0 


- 1.13418526 


0.1784679 


- 1.134185262 


0.1784687866 


- 1.134184(1) 


0.178471(2) 


10.0 


- 0.75694327 


- 0.06154833 


- 0.7569432708 


- 0.0615483234 


- 0.756943(1) 


- 0.061547(1) 


40.0 


- 0.045853 


- 0.0645673 


- 0.045852780 


- 0.0645672604 


- 0.04585286(7) 


- 0.0645673(9) 


400.0 


+ 0.000082 


- 0.002167 


+ 0.00008190 


- 0.0021670 


0.0000818974(3) 


- 0.002167005(3) 



Results for this kinematics on the cut are given in Table 1. The process H 



77 



was investigated before in Ref. |]I6|| . For the master integral under consideration in |16 
all integrations but one could be performed analytically and only the last one had to be 
done numerically (hence the high precision achieved). 

Similarly, high precision is obtained on the cut in our approach, as is demonstrated 
in Table 1. Here, both real and imaginary part of the scalar two-loop H — > 77 integral are 
shown in comparison with the results of Ref. [|16] . We consider the domain q 2 hr < q 2 < 
100q 2 hr , where q 2 hr = Am 2 {m t = 150GeV). For q 2 close to the threshold, the integral has 
a logarithmic singularity, but still we obtain good stability of the approximants, which 
improves to 8-10 decimals up to q 2 = Wq 2 hr and even for q 2 = 100q 2 hr is still excellent. 

For our methods of analytic continuation to work with such high precision, also the 
Taylor coefficients must be known with high accuracy. In general we should know them 
analytically and then approximate them with the desired precision. A good example are 
the coefficients for the H — > 77 decay, which can be represented as rational numbers 
and which for our purpose were approximated with a precision of 45 decimals using 
REDUCE. In general, however, it turns out that the CPU time needed for indices ?i^30 
(in (H) for 1 = m = 0) is of the order of several hours. Moreover in the arbitrary mass 
case the length of the expressions even for lower indices becomes enormous and is getting 
more and more difficult to keep under control. For these reasons it is not possible to 
obtain analytic expressions for the coefficients with a reasonable effort and that is why we 
develop a proper algorithm. In the following three sections we formulate this algorithm 
according to the items specified in Sect. 2. FORM will be a guide for the development 
of the algorithm and a permanent testing tool by comparing results obtained in different 
manners. 
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4. The numerator of the integrand. 

The first step to be done is the formal Taylor expansion of the integral We 
choose the following approach: each scalar propagator with an external momentum p is 
expanded as (p 2 = 0) 

1 = 1 y ( ~2kp V = 1 

(k + p) 2 — mj k 2 — mf \ k 2 — mf J q % ' 

Dealing with the planar and non-planar diagrams (see Fig. 1) simultaneously, as before, 
we can express the Taylor coefficients under consideration in generalization of (5) - (7) 
by 

F n = Df 00n S 1 (k 1 ,p 1 )S2(ki,p2)S 3 (k2,pi)S 4 (k 3 ,p2) \ Pi=0 

n n 

= E ^""V E ^"""V • KAK h, h) . (19) 

u=Q u'=0 

The differential operator -D/oon (0) applied in ( |T9| ) contains two types of operators 
• with rii = n — 2(i — 1) and 

1 — 1 11 LJ 22 

such that the sum of the powers rii + 2(i — 1) = n. In a first step one has to find a formula 
for the application of the n 12 operator. The result depends on n«, the power to which 
this operator is raised, and the partition of scalar products to which it is applied: 

tt(n t ,n l -j 1 ,n l - J 2) = U n £(k lPl )^ {k 2Vl ) n ^ (^ip 2 ) i2 (hp 2 ) n ^ . 
Performing the differentiation with FORM, we find by inspection 



[ m-\-n i 



tt(l,m,n) = J2 t(/,m,n,j)-(A; 1 2 ) / ~( m+ "^(A ;i A; 2 ) m ^(A; 1 A ; 3) n - J (A ;2 A;3) J 

3=0 

with 

t(l,m,n,j) = r f j!(I-m)!(!-n)! ; 



j J \jj ' " l '(l-m-n + j)\ 
Here and in the following we assume inverse powers of factorials of negative ar- 



guments to vanish. Counting the powers of p\ and p 2 in (|T9|), if we wish to calculate 
A™ u i(ki, k2,k 3 ), we need as next 

2(*- 1 )(*-l)!dd((n-z/- Jl ) x kuiv-m + ji) x k 2 ) = □ir 1) (^iPi) n -^ 1 (^i)^ +J1 , 

(20) 

where the function dd of 2(i — 1) arguments is the totally symmetric tensor, contracted 
with n — v — ji vectors k\ and v — rii + ji vectors k 2 ' 

dd{m,\,m-i) = 8(rni,m 2 ) 
dd(mi, m 2 , 7773, 777,4) = S(mi, 777-2)5(7773, 777,4) + S(mi, 1713)6(1712, 777,4) 

+5(7771,7774)5(7772,7773) etc. , 
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where 5 is Kronecker's delta. Similarly for the application of D 22 we write 

2^ 1 \t-l)\dd((n-u' -j 2 ) x k 1 ,(i/-n i + j2) X k 2 ) = {k x p 2 ) n ~ v ' -i\k z p 2 y' ~ n ^* . 

(21) 

The totally symmetric tensor is implemented in FORM 2.2b and can be used for 
performing the differentiation. For high indices n, however, in this manner the differen- 
tiation still requires too much time (for n = 32 appr. 8 hours on the HP735) and the 
expressions are too lengthy for practical use. Therefore the idea is to find an analytic 
expression for (pQ|) and ([H]) and thus to obtain finally a formula for the coefficients a^, 
in (^). Indeed ddim x k±, (I — m) x k 2 ) = dd(m, I — m) can be written as 



dd(m,l-m) = J2 d ( l i m , j)( k i)^ (ki)^ 1 (kik 2 ) j (22) 

j = P(m), 
Aj=2 

with P{m) = (1 - (-l) m )/2 and 

m\(l-m)\ 

d{l,m,j) - 



2 H(t^i)!(2ti)y! ' 

(|22"D has been verified again by inspection of results from FORM's dd- (...) function (in 
FORM 2.2b an algorithm was used but not the evaluation of a formula). 
Summing over all partitions of scalar products, ( |19D yields 

A„ u >(ki, k 2 , k 3 ) 



with 



= <„ + 1) 2- , rt"- 1 ' B £ (-i)'-^r(n + d --i) Si 
^ ' r( n + d - 2)r(n + f ) fct ' m! 2 

min(nj,n— min(nj,n— i/) / \ / \ / /\//\ 

J i=max(O ) n i -i/)ja=max(0,n i -i/') \ Jl / \' h Jl / \ J2 / V tj -' 2 / 

•tt(n i; rij - ji, rij - j 2 ) • dd(n - v - j h v - m + ji) ■ dd(n - v - j 2 , v' - m + j 2 ), 



where the second drf-factor depends on k\ and k 3 instead of k\ and k 2 (see (|19j)). Rewriting 
(|22|) in the following more adequate manner: 



dd(n -v-ji,u-n i +j 1 ) = Y^ 0(cr - [uj - (i - 1)]) 

<T=0 

• d(2(i -l) >n - v -j 1>Vj - 2a) {kiy-^+^klYihhp- 2 - 
with Vj = v — rij + ji we obtain 

M [5] M 

M =0 cr=0 T=0 

.^^-(^O+M^a^^-^-^ 9 (23) 

where the coefficients b^f T are given by 
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= {n + 1) 2 x ~ n - 1 (n-v)\(n- v')\v\ v'\ ^ — ^ j- (24) 



T{n + d-2)T{n + -) 



2 2p L2H 1 J 

S(-4)«(i-i)ir( B + 5 -0 



A' 



1 



min(ni,n— f) min(ni,n— !/') 

V 2 j V 2 

^ \ , ^ / \ (ni — j — p)\(ni — k — p)\(p + j + k — TiAl 

j=max(0,ni-i/) fc=max(0,rtj— i/) \ 1 J " > V 1 r j \r ' J ' »/ 



(a - [j -(„-„) + (<_ 1)])!(^ - 2a)!(r - [* - (n - i/) + (i - 1)])!^ - 2r)! ' 



with A = v + v' — 2fM , p = p — a — t and 
Uj = v — rii + j , v' k — v' — rii + k . 

First of all it is interesting to note that for the planar diagram ( f23|) reduces with 
ks = ki to the simpler form (0) if we set 

£ £ C" ■ (25) 

o-=0 T=0 

i.e. the coefficients a™^,, which in Ref. |l| (appendix A) were only given explicitly as 
rational numbers for a limited number of indices, are now obtained analytically as five- 
fold sums. The nice property of the representation is that, once it has been checked 
against Ref. ]IJ (which has been done with FORM) for the planar diagram, the coefficients 
Ku' aT f° r the nonplanar diagram are checked simultaneously. 



5. Cancellation of the numerator scalar products. 

Since we are here interested in developing an algorithm for the calculation of the 
Taylor coefficients, the above is the first necessary step. The next is to investigate the 
possibility of cancelling the numerator scalar products of integration momenta against 
the bubble propagators Cj (i = 1,...,4). While in a formula manipulating language 
this is done blindly e.g. by using the "repeat" command of FORM, here we have to 
find a detailed prescription if finally our algorithm is to be implemented in terms of a 
FORTRAN program. Cancellation of scalar products yields "genuine" two-loop bubble 
integrals which are investigated in terms of recurrence relations in Sect. 6 and factorized 
one-loop integrals. 

At first we study the planar diagram. Due to (H), © and ([TT| ) we can write 



F„ 



= E 

1 

<2y+v' — 2fi 



n-v n-W 
c l c 2 



C 3 C 4 



v+v'-2n 



(kl + kl-mir^- £ {k\ + kl~ 



m 



6) 



(2k x h 



\a-l 



■ c 6 



a=l,odd 



s 



Let us consider the term (A Q = v + v' — 2// — a; a = 0, • • • , v + v' — 2/i; Ao = A) 

where 5 = a + [3 + 7 was introduced. Accordingly we have to cancel scalar products in 
the following combination 



n-u+l n-u'+l v+lui+1 > 

2 7/i 2 



In both factors of (EBp the sum of powers of k\ and k^, respectively, is larger in the 



denominator than in the numerator so that complete cancellation is possible. Expanding, 
e.g., 



n—11—8 



(*?)"-"-'= E 



- A* - ^ 



K=0 \ K / 

we have to consider the ratio 

V — [l — S — K — 1 
TDK _ <h 

1 .n-u'+l 
c 2 

Depending on the sign ofz^i — 1/ — fi — 8 — k — 1, the partial fraction decomposition 
is performed in the following manner (see also Ref. 0, ch. 10): 

1) v\ > 0: in this case we simply have 

"' f v\ \ (m| — m\) 1 



■ I n—v'—vi+l+i 
i=0 \ b J c 2 



i.e. there remains no c\ in the decomposition. 
2) v x < 0: 

K = E (-i)'' ' 



4 =o v n-i/ 



+ E(-i) kl1 



+ A 1 



8=0 



H - 1 ; (mf - mDM+^r^ 1 ^ 



Similarly we proceed for the fcg-dependent part of (|26|). Expanding 

ft 



K=0 



we deal with 

-iTo 



2 w+1 > 



introducing z/ 2 = /i + 7 — 2/ — ft — 1 and performing the partial fraction decomposition 
as above. Since, however, C5 = k\ — m| is also /^-dependent, for each power of I/C3 and 
I/C4 a further decomposition has to be performed, like e.g. 

* 1 1 

2—1 (^2 ™2\n+l-i J+l ^ 



4 +1 c 5 ~ ^ (m| - m 2 )P+^ 4 +1 (m§ - ml)^ 1 c 5 



9 



For a = we obtain in this manner "genuine" two-loop integrals of the type (^) which 
will be dealt with in Sect. 6, while for a > 1 the factorized one- loop integrals are of the 
type (see also Ref. @) 



(fc? - ra?)" 1 ^ - ml)"* 



with 



<i d /c , , d/ .j, v .^(^ ^ 2X2)4 



2 =^ 7T 2 ( — 171 * ' ' 



Somewhat differently works the cancellation of the numerator scalar products of 
the integration momenta for the non-planar case since here k\ occurs also in C4 = k\ — m\ 
which is raised to high inverse powers. Therefore it is advisable to expand the full product 
of scalar products of different momenta in ( P5| ) in terms of squares of integration momenta, 
i.e. we write with Ai = v — \i — cr + t, A2 = v' — /i + a — r, A3 = fi — a — r and 
A = Ai + A2 = v + v' — 2/i, A = A1 + A2 + A3 = v + 1/ — \i — a — r 

{hk^ihk^ik^ = ^{kl + k\ - k\f\k\ - k\ + k^\k\ - k\ - k 2 3 ) x * 

= qaE(^i) a Poi(k2,k 3 ) . 
Z a=0 

A little algebra also allows to expand P a : 

p a (k 2 2 ,ki) = (-lrEt-i^^y^r^ with (2?) 



min(a, A3) 



fa,P - 


E 




i=max(0,a— A) 


and 






min(k,f3) 


n k - 


E 




£=max(0,k— (a 



9a,p9t-i (28) 



fc^, e.g., can now be completely cancelled only if 

n — fi > A = n — [i + 9 

with = v -\-v f — 72 — (7 — t < 0. Since obviously we do not always have 9 < (in contrary 
to the planar case with n — (u + v') + // > 0, see (|7])) not all fcf can be cancelled, and 
similarly the same holds for k\ and k\. This means that we will obtain integrals of the 
tyP e © with negative indices for which explicit expressions will be given at the end of 
Sect. 6. In fact these integrals are also factorized one- loop integrals which only appear in 
a somewhat different manner than in the planar case. 
The following expansions are needed now: 



n— a— r— a 



K=0 



{klY^ = Ef a + ^C3 +/3 " K (-3) K and 



K=0 

T+a-13 



re=0 
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and partial fraction decompositions must accordingly be performed for the ratios: 

u—a—T—a—K—l a+f3—v—K—l T+a—0—v'—K—l 

- 1 ttt— ,^ and °A . 

c n 2 - u+1 c 5 c 6 

This works out in analogy to the planar case and will not be discussed in further 

details. 

6. Recurrence relations 

As was discussed in Sect. 2, the final step of our algorithm consists in the evaluation 
of the bubble integrals ([J) after the Taylor coefficients are finally expressed in terms of 
these. Their resolution is supposed to be performed in terms of recurrence relations. 
These were first considered in [[17j. One can get such relations from the identity: 

d d hd d k 2 d ( Ak lfl - B k 2il \ _ 
Wl - ™IY dkT, \n - mlYHh - k 2 f - mlY< ) = ° (30) 

with arbitrary constants A and B. Taking, for example, A = (ml + m 2 — m\)/2/m\B we 
obtain a recurrence relation for V#(a,/3,7 + 1) and Vb's with their sum of indices equal 
to a + /3 + 7. For details on recurrence relations for Feynman diagrams see also Refs. [|l8f 



IS]. Explicitly we obtained the following relations for two-loop bubble diagrams with 



three masses (see [20]): 



B 



2(j! - l)mf I 
+ (2(ji-l)+j 3 -d)}V B 



I 2 2 i 2\ 

(m-^ — m 2 + m 3j 



V R 



A(m 1 ,m 2 ,m 3 ) 



{2j 3 m 2 3 2 



B 



32-1 

+ (j 2 -l)(m?-ml-ml)(3--r 

2(j 2 - l)m 2 3 + 2j 3 (ml - mj) + (j 2 - 1 - d){m\ - m\ + ml)] 2~) V B 

A(mi ' m2 ' W3) { 2jl m? l+(3--2-)3- 



(31) 



(32) 



+ 



+ 



J3-1 



+ {j 3 -l){m 2 1 -m 2 2 + ml)(2- -1 



(33) 



2(j 3 - l)m x + 2ji(m 2 - m 3 ) + (j 3 - 1 - d)(m x + m 2 - nr., 

where l ± V B (ji, • • • , m 1 , . . .) = V B (ji±l, . . . , m 1; . . .) etc. and 

1 



A(m l} m 2 ,m 3 ) 



2m\m\ + 2m\m\ + 2m\m\ — m\ — m 2 — m 3 



(34) 



The idea of their application is to reduce all integrals to the master integral 

V B (1, 1, l,mi,m 2 ,m 3 ) (35) 

and some simple tadpole-integrals, which are obtained when one of the indices is zero. By 
inspection one observes that applying all three recurrence relations one after the other, 
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the first index does not increase and the sum of all indices decreases by at least 2 (by 1 
in each of the last two steps). In this manner j 2 or j 3 must get zero and the procedure 
can be stopped. 

A particular point we have to observe is that the integrals (§) may diverge, i.e. writing the 
space-time dimension as d = 4 — 2e, we will have poles in e . Since we intend to develop our 
algorithm in such a manner that it can be implemented in terms of a FORTRAN program 
we have to take special care of these poles, i.e. we have to split the integrals @ into their 
finite and their divergent part. The master integral ( |35D and V B (1, 1, 2, mi, m<2, ma) are 
the only ones which have a pole of second order (for details see also 0), and only the 
integrals V B (1, 1, n, mi, m 2 , m 3 ) with n > 2 have poles of first order, all others being finite. 
Thus we write 

V B (l, 1, n, mi, m 2 , m 3 ) = 1, n, mi, m 2 , m 3 ) + ^ I (1, 1, n, mi, m 2 , m 3 ), n > 2. (36) 



In what follows we will show a path of how to resolve equations (|3T|) - ( p3|). For 
convenience we drop the masses in the argument list. Using recurrence equation (|33|) we 
get 



V B (l,l,n) = 
— {2ml [Vb(2, 1, n - 2) - V B (2, 0, n - 1)] + 



7? 



(n - 1)K - ™2 + [Vb(l,0,n) - Vb(0, l,n)] + 
2(n - l)m^ + 2(m^ - 7713) + (n - 1 - d)(mf + m^ 



Vb(1,1,77 



1)}. 



Here Vb(2, 0, n — 1), Vg(l, 0, n) and Vb(0, 1, n) are "trivial" and will be substituted 
explicitly, e.g. 



V B (0,m,n) 



\m+n 



r(m-|) 



77T 



™/2 r n 



77T 



n/2 



T (ml 



2 

m 2 v 21 



Tin] 



(37) 



with 



r(i 



1 -6. 



(38) 



1 

2/ ~ ~~E 

This indicates the occurrence of a pole in e which gives a contribution to I(l,l,n). 



V B (2, 1, n — 2) is found by an application of recurrence equation (131 

1 



Vb(2,1,77-2) 



2m? 



{(77-2) [Vb(0, 1, 77 - 1) - Vb(1, 0, 77 - 1)- 



(m?-mi + m§)Vb(l,l,n-l)] + (77 - c?)Vb(1, l, 77 - 2)} 



Therefore we can inductively calculate Ve(l, l,n) for all 77 > 3 using the master integral 
V B (1, l, l). Care has to be taken in evaluating terms of the form d V B (1, l, 77 — l) on the 
r.h.s. which yield a finite and an infinite part with d = 4 — 2e . Using recurrence equation 
(|32"D, we get moreover 

V B (l,m,n) = 

12777773= [Vb(1, m - 2, 77 + l) - Vb(0, m - l, 77 + l)] + 
m — 1 1 

(m — l)(m^ — wij — 7773) [Vb(1, m, 77 — l) — V B (0, m, 77)] + 
2(m — l)m 3 + 2n(m\ — ml) + (m — l — d)[m\ — m\ + 777,3) Vb(1, m — l,n)| 
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so that V B (l,m, n) is computable, and finally an application of recurrence equation (|3"T|) 
leads to 



V B (l,m,n) = 

— - {n [V B (l - 2, m, n + 1) - V B (l - 1, m - 1, n + 1)- 




(m? - + m|)F B (i - 1, m, n + 1) + (2(1 - 1) + n - d)V B (l 




Implementing the recurrence relations numerically, it is advisable, in order to take 
properly into account the poles in e, to produce explicit relations for the lower indices 
by means of FORM. For the V B (1, 1, n), e.g., the Jjj pole from the master integral causes 
"extra" divergences and for the V B (l, m, n) divergences come in for lower indices whenever 
a V B (1, 1, n) is encountered in a recurrence relation while for higher indices these integrals 
are finite as can also directly be seen from (|9|). 

Apart from that, the recurrence relations are implemented in exactly the order as 
described above: first V B (1, 1, n) for all n, secondly V B (1, m, n) and finally V B (l, m, n). In 
each of these cases the recurrence has been implemented for all I + m + n < J (I, m, n > 1; 
for negative indices see below), where the large integer J is to be chosen according to the 
number of Taylor coefficients needed. 

Of course one may have doubts about the numerical stability of such a recursive 
approach, in particular if one knows that indeed the Taylor coefficients must be calculated 
with high precision as mentioned in Sect. 2. Ordinary FORTRAN double precision will 
clearly not do the job. We used the multiple precision package written by D.H.Bayley 
which also provides an automatic translator for any FORTRAN program. The 
requested precision is here defined at the beginning of the program. With this package 
we used for J = 62 and 100 decimals precision 42 seconds on a HP735 (49 seconds for 150 
decimals precision). The numerical results were tested against numerics performed with 
REDUCE for explicit expressions of the integrals in the equal mass case. In principle, 
however, the best precision test will always be to increase the number of decimals used 
in the calculation. In this manner we can be sure to have an effective algorithm for the 
calculation of two-loop integrals. 

To conclude this sections we give an explicit formula for V B (a, /3, 7) for the case 
that one of the indices is negative. In this case the V B s can be reduced to factorized one- 
loop integrals, which makes the application of the above recurrence relations superfluous. 



V B (a,(3, 1 ) = (-lY a+ ^ 



I 



d d h d d k 2 pi - k 2 ) 2 - 
(k\ - m 2 l ) a (kl - m. 




.2117 



(39) 



M) 2 |7|! E 



[|7l/2] (d 



a-l /3-1 



T(l+r-l-l)T(l + q-l-l) 
r\q\(a-l - r)\(/3 - 1 - q)\ 



l\ 



E E 



r=max(0,2l+a— 1) 9=0 




(| T | -2l-a- P + r + q + 2)\ 



where 7 < and s = —ml — m\ + m\. 
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7. Conclusions 



An effective method to calculate Feynman diagrams has been developed in [p]]. 
In the present work we have been able to work out details of an algorithm, which will 
finally allow to elaborate the above method into a "package" for the evaluation of two- 
loop three-point functions and possibly beyond. The essential new points worked out here 
are described in Sects. 4 and 6: deriving an explicit formula for the numerators in the 
integral representation of the Feynman diagrams's Taylor coefficients and demonstrating 
the possibility to evaluate the recursion relations for the bubble integrals numerically by 
means of the multiple precision package of [^TJ. In the development of this algorithm 



FORM has been the main tool in the formulae manipulation. 
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